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' We study numerically evolution of metric perturbations during reheating in 

q_) ■ a model with two fields and a strong parametric resonance. Our calculation 

1^ . is fully nonlinear and includes gravity but is restricted to spherical symmetry. 

OO \ In this model, super-Hubble metric perturbations can grow during reheating 

\ only due to effects nonlinear in fluctuations of the fields. We find that they 

C^ ' indeed grow and, soon after the growth begins, dominate variances of the 



> 



o 



metric functions. Thus, the metric functions become smooth but varying 



Q> . significantly over large scales. Their profiles at late times are interpreted as 



signalling a gravitational instability and formation of a black hole. 
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Inflationary cosmology (reviewed in books |]T|) provides a very attractive explanation of 
the origin of the large-scale structure of the universe 0. If the only important interaction 
of the scalar fleld (the inflaton) is gravity, metric perturbations at super-Hubble scales can 
grow during inflation and during the transition to radiation-dominated (RD) era, but not 
during the RD era itself. Thus, they are preserved until the time when a gravitational 
instability is ready to convert them into the observed structure. 

The question remains, however, whether that scenario can be changed by additional 
interactions of the inflaton (a self- interaction or interactions with other flelds). Speciflcally, 
right after the inflation, at reheating, the inflaton still has a large homogeneous component, 
and it is easy to imagine that fluctuations of various flelds (and essentially of any scale) 
produced at that time will kick particles out of the homogeneous modes, into modes with 
large, super-Hubble wavelengths. 

A natural arena to address the above question is inflationary models in which reheating 
starts out explosively, via a parametric resonance 0. In rigid spacetimes, the growth of large- 
wavelength modes due to rescattering |^, and the resulting formation of smooth (turbulent) 
power spectra, are clearly seen in a variety of inflationary models ^j^- To extend these 
analyses to super-Hubble scales, it is necessary to include gravity, i.e. perturbations of the 
metric. 

The influence of parametric resonance on the evolution of metric perturbations at flrst 
order in fleld fluctuations was addressed in 0-0 . While it is possible to identify models in 
which the growth of super-Hubble metric perturbations occurs at that order P,0, in many 
popular models such growth (if any) can only appear at the second order (e.g. through the 
scattering process described above). Such second-order effects are quite generic. Fortunately, 
models with a parametric resonance are amenable to a fully nonlinear numerical study, based 
on the classical approximation [Q. 

Here we present results of such a study, including gravity but restricted to spherical 
symmetry, for the model of two interacting fleld with the potential 

\/(0,x)=/0V/2 + A0V4, (1) 

where (j) is the inflaton, and x is a fleld that initially has no homogeneous component. Both 
flelds are minimally coupled to gravity. 

We choose one the "magic numbers" for g'^/X, for which the resonance is especially strong 
\T^■. g^/X = 8. With this choice, the resonance peak is at fc = 0, so super-Hubble fluctuations 



of X will start growing right away. Because x initially has no homogeneous component, that 
growth will induce fluctuations of the metric that are still of the second order in 6x- We do 
not believe the choice of g^/X to be of great importance (although a weaker resonance or 
a smaller initial Sx will delay the onset of nonlinearity in fleld fluctuations); we tried other 
values, with very similar results ||12[. Note that the scattering mechanism described above 
will work even in the absence of any preexisting fluctuations of x at super-Hubble scales. 

We flnd that (i) super-Hubble perturbations of metric indeed grow, and (ii) proflles of the 
metric functions eventually undergo a signiflcant change, which we interpret as development 
of gravitational instability and formation of a black hole. 

The flrst of these results is in agreement with a result obtained by Easther and Parry 
for the X(f)^ model in planar geometry [|T^, but the second one is not. We do not think 
that the difference is due to the different choice of a model, but are inclined to attribute 



it either to the difference in geometry or to the presence of an additional radiation field in 
their calculation. 

When we speak about terms of the second order in 5xi we envision fluctuations of the 
energy density of the form 

5E^ oc j d^pF{p, k)5xp<5Xp+k ' (2) 

with nonzero k {F is some function; a corresponding term with k = does not contribute 
to 6E, but only to the average energy). Terms like (Q) appear both directly from (^x)^ 
terms in the energy, and from terms proportional to 50, after a sizable S(j) is generated by 
rescattering. 

We caution, however, against relying on the size of 5E as the criterion for formation of 



black holes, as seems to be the trend in recent literature IjT^. One needs to solve explicitly 
for the metric functions. In Newtonian gravity, the gravitational potential will typically be 
smoother that the mass distribution. Similarly, in our simulations (using general relativity), 
we find that while the profiles of the fields and x (and of the spatial curvature), at the time 
when a black hole forms, are fairly rugged, those of the metric functions are exceptionally 
smooth. We stress that at that time the power spectra of the fields and of the metric 
functions are determined by rescattering, rather than by their initial spectra or the original 
shape of the resonance. Thus, rescattering becomes a crucial effect. 
We use a spherically symmetric metric in isotropic coordinates: 

els'" = -N^it, r)dt^ + ^\t, r) [dr'' + r't/fi'] . (3) 

Equations of motion are written in the Hamiltonian form, the Hamiltonian pairs being 
(0, 7r<^), (x, vr^), and (-0, K), where K is the trace of the extrinsic curvature: K = —Gip/Nip. 
The lapse N does not have a pair; it satisfies an ordinary differential equation that involves 
only spatial derivatives. The energy and momentum of the scalar fields are 

E = inl + 4)/2^i^ + (0'^ + x")/2^' + V{<P, x) , (4) 

P=(7r^0' + 7r^x')/^', (5) 

and the energy and momentum constraints are 

Ce = R + SV't/'/V'^ = , (6) 

Cm = P- 2K'/3k^ = , (7) 

where R is the curvature of the spatial sections: 

R/2 = K^E - K^S , (8) 

and K^ = SttG. 

We take the equation for K in the form 



Adding the constraint function Ce with some coefficient w > supphes the equation with 
a second spatial derivative of ip and, for a range of w, prevents a short-scale numerical 
instability. We present results for w = 3; results for w = 1 are similar. All other equations 
are standard, and we do not write them here. 

It is convenient to work in rescaled variables, in which and x ^^^e measured in units 
of 0(0), the field at the end of inflation, while time and distance are measured in units of 
1/\/A(/)(0). Comoving momenta k are in units of vA0(O). The evolution equations retain 
their form, with the following replacements: A ^ 1, ^f^ — ^ 5'^/A, and k^ — > fi;^0^(O) = 3.08. 
The numerical value for 0(0) is adopted from the previous work ^. 

In the classical approximation, fields are regarded as classical entities with random initial 
conditions. Their Fourier components in the initial state are Gaussian and distributed 
according to certain power spectra. We choose these initial spectra as follows. 

Amplification of quantum fluctuations of the metric during inflation results in an approx- 
imately scale-invariant spectrum for if): (dV'kl^)) "-^ 1/^^; double brackets denote averaging 
over the quantum state. Then, from the energy constraint, ((|7r,|^p)) oc A; at large k. We 
choose to cut off the spectra at large k to remove the corresponding large contribution to 
energy, so it does not cause spurious gravitational effects at early times. Thus, we use for -0 
an initial power spectrum of the form k~^ exp(— /c^/ju,^), with /i^ = 10. This cutoff is imposed 
only at the initial instant, and the results at large times are insensitive to it. The overall 
magnitude of the -0 spectrum is chosen to get a realistic size of fluctuations on super-Hubble 
scales, (5'0/'0 ~ 10~^. 

Initial spectra of 0, Xi ^'^^ ^ cire then obtained as follows. Fluctuations of x ^^e placed 
in adiabatic vacuum, with a cutoff factor as above. This leads to an initial ^x ~ 3 x 10~^Mpi. 
Fluctuations of are determined from the constraints. We either solve the constraints to 
the first order, or use a trivial exact solution: 50 = 0, with 5K determined from (|^), and 5ti^ 
from d^). Results for these two types of initial conditions are quite similar; here we present 
those for the exact solution. 

The Hamiltonian structure of the equations suggests that it is suitable to use a leapfrog 
algorithm, in which the coordinates are viewed as defined at full time steps, and the momenta 
at half-steps. The lapse N is updated from its differential equation at full steps and is 
extrapolated to half-steps. The equation for A^ does not fix its overall normalization, so A^ 
can be arbitrarily scaled at each spatial section to define a convenient time variable. We use 
conformal time, obtained by setting A^ at the outer end of the grid equal to the scale factor 
(defined below). We use the boundary condition /„ = /„_i at the outer end; / is any field 
or metric function. 

At times shown in the figures, -0^, the square of -0 at the end of the grid, grows steadily 
with time and even at the latest times deviates relatively little from a linear t dependence 
corresponding to an RD universe. Motivated by that, we adopted the following definitions 
of the scale factor a{t) and the comoving Hubble parameter H{t): a = -0^, H = — aA'„/3. 

We output two types of quantities as functions of the conformal time. Some are built 
from ordinary lattice averages, such as 



n-l 



if) = s^' E ^V. , (10) 



i=l 



where 5*2 = {2n'^ —3n'^ +n) /Q. These include ordinary root mean square (r.m.s.) deviations- 



square roots of the variances (((5/)^), where Sf = f — {f). The other type of quantity is 
r.m.s "superdeviations" cr(/), which measure the super-Hubble content of a field: 

a\f) =471 f dkPPfik) , (11) 

•^ "'mill 

where Pf is the power spectrum oi 6f = f — (/), and fcmin = tt/L; L is the radius of the box. 

Here we present results for L = SOOvr and n = 8192. In Fig. 1 we plot two quantities 
characterizing the growth of super-Hubble perturbations: a{R)/2n'^{E), sud a {ip) / {ip). Both 
clearly grow, and that growth begins around the same time as the growth of 5(j). This is of 
no surprise since all these effects are of the same, second order in 5x- (The resonant growth 
of 5x begins much earlier, at t ~ 10, and ends at t ~ 80.) We also plot the divided by ('0) 
ordinary r.m.s. deviation of ip, which includes all modes. We see that at t > 60 the full 
r.m.s. deviation is barely distinguishable from (j{il))., i.e. at these times almost all of the 
variance of ip is in super-Hubble modes. Same is true for A^, but not for example for R. 
This quantifies our conclusion that at times of interest the metric functions are smooth. 

To judge the quality of these results, we plot in Fig. 2 the superdeviations of the 
constraint functions Ce and Cm, relative to cr(i?) and cr(P), respectively. The large spikes 
in a{CM)/<y{P) occur at times when the homogeneous component of cf) reaches a maximum 
or minimum, and cr(P), by which we divide, becomes especially small. 

In Fig. 3 we plot profiles of the lapse function N{r) at different moments of time. The 
values of t give a good approximation to 1/ H, the comoving Hubble radii at these times 
(with H defined as above), and one sees that the variation of N mostly occurs over super- 
Hubble scales. We interpret A^ crossing zero as formation of a horizon of a black hole (or, 
more precisely, what will look as a black hole to an observer at large r): a point where 
A^ = 0, and if) is finite, cannot be reached by light from any other point in a finite time. 
For a Schwarzschild black hole in asymptotically flat space, the lapse function, in isotropic 
coordinates (^), crosses zero at r = ro/4, where tq is the Schwarzschild radius. In that 
case, the isotropic coordinates cover only the region outside the horizon of the black hole, 
but they cover it twice. In our calculation, the spacetime is not asymptotically flat, and 
at present we are unsure how to interpret the "inside" region, where iV < (see, however, 
remarks below). 

We find the presence of super-Hubble modes essential: if we reduce the size of the 
integration box, so that no such modes are left by the time when nonlinearity becomes 
important (by choosing e.g. L = AOn), N does not cross zero. For a slightly larger box 
{L = TOtt) we observe formation and subsequent disappearance ("evaporation") of a black 
hole. These results correspond well with our observation that Sip/ip reaches smaller values 
in smaller boxes. 

In a large volume, after the zero of A^ first occurs, it propagates to larger r, until it 
becomes comparable to the size of the box, and the constraints break down. Shortly before 
that, however, we observe an interesting change in the profiles of the fields, see Fig. 4. The 
field in the "interior" region becomes smooth and reminiscent of an excited state of a 
soliton star, cf. ref. [0, while the field x is pushed out of the "interior" (we suppose, by its 
repulsive interaction with 0). 

In the full three dimensions (as opposed to the spherically symmetric case) we expect 
formation of many black holes in different regions of space. In addition, gravitational waves 



(absent in spherical symmetry) may be produced in these inhomogeneous spacetimes even 
more efficiently than indicated by calculations in rigid spacetimes |T^. It remains to see 



how these various effects constrain inffationary models in which they occur. 
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FIG. 1. Growth of super-Hubble perturbations: normalized superdeviations (see text for 
the definition) of the metric function ip (sohd hne) and the spatial curvature R (long dashes) as 
functions of the conformal time. The normalized full r.m.s. deviation of ip, including all modes, is 
shown for comparison (short dashes). 
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FIG. 2. Normalized superdeviations of the constraint functions Ce (solid line) and Cm (dashed 
line) as functions of conformal time. 
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FIG. 3. Profiles of the lapse function N at several moments of time. We interpret the crossing 
of zero as formation of a black hole. 
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FIG. 4. Profiles of the fields (p (solid line) and x (dashed line) at t = 85. 



